Bare Necessities

Author

Adam Kemberling

Published

December 3, 2024

Code
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
library(nlme)
library(ggimage)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Everything in this markdown uses the size spectra exponent values where I shifted the wmin parameter to 16g. The wmin values are fixed at 16g across all groups.

Code
# # Data used for Wigley estimates
# wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>%   
#   left_join(area_df)
# 
#### 1g Cutoff Data  ####
# # Load the three regional datasets
# wigley_bmspectra <- read_csv(
#   here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
#   mutate(survey_area = factor(survey_area, levels = area_levels))
# 
# 
# # Load Shelfwide
# wigley_bmspectra_shelf <- read_csv(
#   here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
# 
# 
# 
# # Join them together
# wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
#   mutate(metric = "bodymass_spectra",
#          community = "Wigley Subset",
#          survey_area = factor(survey_area, levels = area_levels_all)) %>% 
#   left_join(area_df) %>% 
#   mutate(area = factor(area, area_levels_long_all))


# Modeling datasets

# 
# # # Model Datasets for Wigley Species Subset
# wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))
# 
# # load shelfwide stuff too
# wigley_bmspectra_shelf <- read_csv(
#   here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
# 
# # Bottom Temperatures
# bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))
# 
# 
# # Bolt that on
# wigley_bmspectra_model_df <- bind_rows(
#   wigley_bmspectra_model_df,
#   left_join(
#     wigley_bmspectra_shelf, 
#     bot_temps, 
#     by  = join_by("est_year" == "year", survey_area, season))) %>% 
#   left_join(area_df) %>% 
#   mutate(area = if_else(survey_area == "Northeast Shelf", "Northeast Shelf", area),
#          area = factor(area, levels = area_levels_long_all),
#   season = fct_rev(season))
Code
# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>%   
  left_join(area_df) 

# Data used for 16g wmin
wigley_in_16 <- wigley_in %>% 
  filter(wmin_g >= 16)

# Load the 16g minimum size cutoff
wigley_bmspectra_model_df <- read_csv(
  here::here("Data/model_ready/wigley_community_16g_bmspectra_mod.csv")) %>% 
  mutate(
    metric = "bodymass_spectra",
    community = "Wigley Subset",
    survey_area = factor(survey_area, levels = area_levels_all),
    area = factor(area, area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))
  )

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

Temperature and Size Structure

There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

  • Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
  • James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

The metabolic theory of ecology provides a basis for “why” the above rules might be observed.

  • Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Seasonal Movements and Metabolic Efficiency: MTE details how metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

This theory would suggest that larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and to reduce the energetic costs of thermoregulation in warmer periods.

From what I see at first-glance the size distribution along the shelf is responding consistently with temperature changes (Spring to Fall) on shorter time scales than hypotheses that predict temperature related changes in growth.

Code
wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  ggforce::geom_mark_ellipse(aes(fill = season), alpha = 0.1) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(
    shape = guide_legend(nrow = 2),
    fill = guide_legend(nrow = 2),
    color = guide_legend(nrow = 2)) +
  theme(legend.box = "horizontal") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "The broad relationship between temperature and size distribution is bolstered by\nseasonal community changes in alignment with seasonal temperature differences.",
    shape = "Survey Area",
    fill = "Season",
    color = "Season",
    x = "Bottom Temperature")

Importance of Seasonality

The extent that the regional spectra trends follow a negative linear relationship to increasing temperature is largely due to seasonal differences.

We know that individuals may migrate and/or occupy different habitats seasonally, leading to situations where the community compositions from spring to fall are an apples to oranges comparison. Seasonal movements may be be partially explained by expectations from the MTE, but they also could be separate ecological processes like seasonal foraging or spawning opportunities.

Either way, seasonal differences are large and they justify treating survey seasons as separate groups with independent trends.

Anomalies Model

I wanted to try and avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

  • Season
  • Survey area
  • Bottom temperature anomaly

Error structure: AR1 autocorrelative errors

By using the season * region interaction groups, we can characterize the typical bottom temperature for each combinations using the long-term seasonal average, and quantify warming using departures from those long-term averages.

The model looks like this as an equation (90% sure), where \(\alpha\) is the intercept, \(\beta x_{tjs}\) is our coefficient for a continuous independent variable (ex. bottom temperature), and \(e_{tjs\)` is our error:

\[y_{tjs} = \alpha + \beta x_{tjs} + e_{tjs}\] and our error depends on the previous timestep:

\[e_{tjs} = \phi e_{t-1,js} + w_{tjs}\] With \(w_{tjs}\) adding noise.

\[w_{tjs} \sim N(0, \sigma^2) \] All of this happening for each \(year_t\) within \(region_j\) & \(season_s\).

I need to consult the Zuur paper about equations in ecology again, but this seems close.

Code
# Just drop the random effect part for now
anom_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(anom_model_ar)$corStruct


# # check the model
# check_model(anom_model_ar)
# plot(check_collinearity(anom_model_ar))
# plot(
#   check_collinearity(
#     nlme::gls(
#       b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), 
#       data = wigley_bmspectra_model_df, 
#       correlation = corAR1(form = ~ est_year | area/season)
# ))
#   )


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    anom_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = anom_model_ar,
  specs =  ~ area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop") %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in:
1. GOM-Spring
2. GB-Spring
3. GB-Fall
4. SNE-Fall

Bottom Temperature Model

We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be than the anomaly model.

I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.

Code
# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct
# 
# # check the model
# check_model(temp_model_ar)

# # Collinearity without interactions
# plot(
#   check_collinearity(
#      nlme::gls(
#         b ~ area + season + scale(bot_temp), 
#         data = wigley_bmspectra_model_df, 
#         correlation = corAR1(form = ~ est_year | area/season))
# ))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ area * season,
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp_results <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Seasonal Bottom Temperature")
local_btemp_results

Verdict: Temperature itself performs better, and is simpler to produce/explain.

Code
data.frame(
  "Model" = c("Temperature Model", "Anomalies Model"),
  "AIC" = c(AIC(temp_model_ar), AIC(anom_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()
Model AIC Delta-AIC
Temperature Model -283.6706 0.00000
Anomalies Model -261.7763 21.89431

Thoughts:

Consistent seasonal shifts in the Shelf-wide spectra suggests to me that the community is more responsive to within-year time scale factors and maybe less sensistive to long-term trends in temperature than I originally expected. This also raises questions about how the spring/fall size distribution differences are being achieved and whether to view them as expressions of the same “community” or separately.

From the decadal variability work we have reason to believe that individual species are both not (perfectly) tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

In light of that knowledge, there must be some compensatory mechanisms happening across this broader community for the temperature+spectra relationship to be so limited/situational.

Region * Season Pairwise Posthoc

The following figures and tables compare marginal means and post-hoc comparisons for the model above using bottom temperature.

What are the differences across spatial scales and seasonally?

Code
# Region and seasonal posthocs
regseas_phoc <- emmeans(
  temp_model_ar, 
  list(pairwise ~ area + season), 
  adjust = "tukey",
  type = "response")



# Custom Plot
(intercepts_emmeans <- regseas_phoc$`emmeans of area, season` %>% 
    as_tibble() %>%
    ggplot() +
    geom_pointrange(
      aes(area, emmean, ymin = lower.CL, ymax = upper.CL, color = season),
      position = position_dodge(width = 0.25), alpha = 1) +
    #coord_flip() +
    scale_color_gmri() +
    labs(
      y = "Marginal Mean: Mass Spectra Slope",
      x = "Area",
      title = "Wigley Species",
      color = "Season",
      subtitle = "Group marginal means plot"))

Code
# Table for pairwise test results
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  arrange(p.value) %>% 
  gt() %>% 
  tab_header(title = "Factor Groups Pairwise Post-Hoc") %>% 
  tab_style(
    style = list(
      cell_fill(color = "#9AFF9A"), 
      cell_text(color = "black")),
    locations = cells_body(
      rows = p.value <= 0.05) 
  )
Factor Groups Pairwise Post-Hoc
Contrast estimate SE df t.ratio p.value
Northeast Shelf Spring - (Mid-Atlantic Bight Fall) 0.943 0.130 471.630 7.252 0.000
Gulf of Maine Spring - (Mid-Atlantic Bight Fall) 0.878 0.134 475.331 6.549 0.000
Georges Bank Spring - (Mid-Atlantic Bight Fall) 0.862 0.150 474.065 5.731 0.000
Southern New England Spring - (Mid-Atlantic Bight Fall) 0.931 0.119 474.694 7.805 0.000
(Mid-Atlantic Bight Spring) - (Mid-Atlantic Bight Fall) 0.978 0.109 473.235 9.014 0.000
Northeast Shelf Fall - Georges Bank Fall -0.384 0.082 474.790 -4.660 0.000
Northeast Shelf Fall - (Mid-Atlantic Bight Fall) 0.733 0.117 439.871 6.287 0.000
Gulf of Maine Fall - (Mid-Atlantic Bight Fall) 1.025 0.118 442.826 8.695 0.000
Georges Bank Fall - (Mid-Atlantic Bight Fall) 1.117 0.122 477.971 9.173 0.000
Southern New England Fall - (Mid-Atlantic Bight Fall) 0.742 0.125 476.724 5.913 0.000
(Mid-Atlantic Bight Spring) - Northeast Shelf Fall 0.245 0.061 458.918 4.011 0.003
Georges Bank Fall - Southern New England Fall 0.375 0.094 475.495 3.976 0.003
Northeast Shelf Fall - Gulf of Maine Fall -0.292 0.077 460.546 -3.815 0.006
Gulf of Maine Fall - Southern New England Fall 0.283 0.089 471.001 3.173 0.051
(Mid-Atlantic Bight Spring) - Southern New England Fall 0.237 0.077 471.559 3.091 0.065
Southern New England Spring - Northeast Shelf Fall 0.198 0.079 470.459 2.518 0.262
Gulf of Maine Spring - Georges Bank Fall -0.239 0.106 477.157 -2.261 0.417
Northeast Shelf Spring - Northeast Shelf Fall 0.209 0.094 463.762 2.226 0.441
Southern New England Spring - Georges Bank Fall -0.186 0.086 477.144 -2.156 0.489
Southern New England Spring - Southern New England Fall 0.190 0.091 474.120 2.079 0.543
Georges Bank Spring - Georges Bank Fall -0.255 0.126 476.856 -2.031 0.578
(Mid-Atlantic Bight Spring) - Georges Bank Fall -0.139 0.070 476.051 -1.970 0.621
Northeast Shelf Spring - Southern New England Fall 0.201 0.105 469.083 1.917 0.657
Northeast Shelf Spring - Georges Bank Fall -0.174 0.100 471.737 -1.739 0.773
Gulf of Maine Spring - Northeast Shelf Fall 0.145 0.100 473.576 1.454 0.909
Gulf of Maine Spring - Gulf of Maine Fall -0.147 0.101 476.177 -1.453 0.909
Georges Bank Spring - Gulf of Maine Fall -0.163 0.122 476.375 -1.340 0.944
Gulf of Maine Spring - Southern New England Fall 0.136 0.110 471.249 1.243 0.965
Southern New England Spring - Gulf of Maine Fall -0.094 0.081 477.748 -1.165 0.977
Gulf of Maine Spring - (Mid-Atlantic Bight Spring) -0.100 0.090 473.850 -1.112 0.983
Gulf of Maine Fall - Georges Bank Fall -0.092 0.084 444.940 -1.092 0.985
Georges Bank Spring - Northeast Shelf Fall 0.129 0.121 469.365 1.066 0.988
Georges Bank Spring - (Mid-Atlantic Bight Spring) -0.116 0.113 477.514 -1.032 0.990
Georges Bank Spring - Southern New England Fall 0.120 0.129 431.576 0.930 0.995
Northeast Shelf Spring - Gulf of Maine Fall -0.083 0.096 469.882 -0.864 0.997
Southern New England Spring - (Mid-Atlantic Bight Spring) -0.047 0.066 477.952 -0.711 0.999
(Mid-Atlantic Bight Spring) - Gulf of Maine Fall -0.047 0.063 470.330 -0.738 0.999
Northeast Shelf Spring - Gulf of Maine Spring 0.064 0.115 477.403 0.560 1.000
Northeast Shelf Spring - Georges Bank Spring 0.081 0.134 475.084 0.604 1.000
Northeast Shelf Spring - Southern New England Spring 0.011 0.097 408.990 0.115 1.000
Northeast Shelf Spring - (Mid-Atlantic Bight Spring) -0.036 0.084 474.661 -0.427 1.000
Gulf of Maine Spring - Georges Bank Spring 0.016 0.138 433.350 0.119 1.000
Gulf of Maine Spring - Southern New England Spring -0.053 0.103 477.368 -0.517 1.000
Georges Bank Spring - Southern New England Spring -0.070 0.123 476.681 -0.564 1.000
Northeast Shelf Fall - Southern New England Fall -0.008 0.088 417.718 -0.096 1.000
Code
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  filter(p.value < 0.05)  %>% 
  separate(Contrast, into = c("group1", "group2"), sep = " - ") %>% 
  mutate(
    group1 = str_remove_all(group1, "[(]|[)]"),
    group2 = str_remove_all(group2, "[(]|[)]"),
    statement = 
      if_else(estimate > 0, " is shallower than ", " is steeper than ")) %>% 
  select(group1, statement, group2) %>% 
  arrange(statement, group1) %>% 
  gt()
group1 statement group2
Georges Bank Fall is shallower than Southern New England Fall
Georges Bank Fall is shallower than Mid-Atlantic Bight Fall
Georges Bank Spring is shallower than Mid-Atlantic Bight Fall
Gulf of Maine Fall is shallower than Mid-Atlantic Bight Fall
Gulf of Maine Spring is shallower than Mid-Atlantic Bight Fall
Mid-Atlantic Bight Spring is shallower than Northeast Shelf Fall
Mid-Atlantic Bight Spring is shallower than Mid-Atlantic Bight Fall
Northeast Shelf Fall is shallower than Mid-Atlantic Bight Fall
Northeast Shelf Spring is shallower than Mid-Atlantic Bight Fall
Southern New England Fall is shallower than Mid-Atlantic Bight Fall
Southern New England Spring is shallower than Mid-Atlantic Bight Fall
Northeast Shelf Fall is steeper than Gulf of Maine Fall
Northeast Shelf Fall is steeper than Georges Bank Fall

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

Code
# Plot annual totals
wigley_bmspectra_model_df %>% 
  filter(season == "Spring") %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.8) +
  geom_ma(
    aes(color = "5-Year Moving Average"), 
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries", linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

Code
# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ area * season + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(est_year > 1978),
  correlation = corAR1(form = ~ est_year | area/season)
)



# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Code
# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

This is our old friend “higher landings when spectra was shallower” that we’ve seen before.

Temperature and Landings Together

If we want to just throw them both in because we’ve done what we can and want to hold them both in the frame at once we get this:

Code
# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  # Model
  model = b ~ area * season * scale(bot_temp) + area * scale(log(total_weight_lb)), 
  
  # Data
  data = wigley_bmspectra_model_df %>%
    drop_na(total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  
  # Autocorrelation
  correlation = corAR1(form = ~ est_year | area/season))

# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp) +  scale(log_land),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%
#                  filter(area != "Northeast Shelf") %>% 
#           mutate(log_land = log(total_weight_lb)),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct


# temp model is much better
# AIC(temp_model_ar)
# AIC(full_model_ar)

Before proceeding it is worth noting that this model is less parsimonious than a temperature or anomaly only model:

Code
data.frame(
  "Model" = c(
    "Temperature Model", 
    "Anomalies Model", 
    "Landings Model", 
    "Landings and Temperature"), 
  "AIC" = c(
    AIC(temp_model_ar), 
    AIC(anom_model_ar),
    AIC(f_model_ar),
    AIC(full_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()
Model AIC Delta-AIC
Temperature Model -283.6706 0.00000
Anomalies Model -261.7763 21.89431
Landings Model -230.3029 53.36772
Landings and Temperature -159.7297 123.94097
Code
# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area!= "Northeast Shelf"),
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp

Code
# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

Code
zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

Cold Pool Index

There are several cold pool indices available that might be helpful in adding context to changes in MAB. Below are the Dupontavice et al. 2022 methods that are used in the SOE report.

The Cold Pool Index (Model_CPI) was adapted from Miller et al. (2016). Residual temperature was calculated in each grid cell, \(i\), in the Cold Pool domain as the difference between the average bottom temperature at the year \(y\) (\(T_{i,y}\)) and the average bottom temperature over the period 1972–2019 (\(\bar{T}_{i, 1972-2019}\)⁠⁠) between June and September. \(CPI\) was calculated as the mean residual temperature over the Cold Pool domain such that

\[CPI_y~=~\frac{\sum_{i=1}^{n}~(T_{i,y}-\bar{T}_{i,~1972-2019})}{n}\]

where n is the number of grid cells over the Cold Pool domain.

Takeaway:

\(CPI\) was calculated as the mean residual temperature over the Cold Pool domain

From this index we can see that the temperatures over the cold pool domain (which include the MAB) have increased over time.

Code
cpool <- ecodata::cold_pool
cpool %>% 
  filter(!str_detect(Var, "se_")) %>% 
  ggplot(aes(Time, Value)) +
  geom_line() +
  facet_grid(Var~EPU, scales = "free")

Code
# # There's an SF
# cpool_sf <- ecodata::cold_pool_sf %>% 
#   mutate(id = row_number())
# 
# ggplot(cpool_sf) + geom_sf(aes(fill = id), alpha = 0.2)

Driver Correlations